Interlayer screening effect in graphene multilayers with ABA and ABC stacking 
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We study the effect of perpendicular electric fields on the band structures of ABA and ABC 
graphene multilayers, and find that the electronic screening effect is significantly different between 
them. In ABA multilayers, the field produces a band overlap and gives a linear screening, while 
in ABC multilayers, in contrast, it opens an energy gap in the surface-state band at low energy, 
leading to a strong screening effect essentially non-linear to the field amplitude. The energy gap of 
a large ABC stack sharply rises when the external field exceeds a certain critical value. 



I. INTRODUCTION 

Recent experimental realizations of atomically-thin 
graphene systems*— open up possibilities of exploring 
their exotic electronic properties. In multilayer films 
composed of more than two graphene layers, the inter- 
layer coupling strongly modifies the linear dispersion of 
monolayer graphene, resulting in various electronic struc- 
tures depending on the number of layers, N£~— The 
band structure can also be changed by applying a gate 
electric field perpendicular to the layer, through gen- 
erating an interlayer potential asymmetry In bilayer 
graphene, for example, an energy gap opens between the 
conduction and valence bands in presence of gate elec- 
tric field£ ' 10 i 12 ' 20 i 23 ~ — and it was actually observed in 
transpor t 28 ' 29 and spectroscopic measurement s 5 ' 6 ' 30 " —. 

In nature, there are two known forms of bulk graphite 
called ABA (AB, hexagonal, or Bernal) and ABC (rhom- 
bohedral) with different stacking manners as shown in 
Fig. [TJ The ABA phase is thermodynamically stable and 
common, while it is known that some portion of natu- 
ral graphite takes the ABC form^ For ABA graphite, 
the effective mass model describing the electronic prop- 
erty was developed for the bulk system 3 ^—, and also 
for few-layer systems.— ~ 10 ' 12 ~ —. The energy dispersion of 
the multilayer graphenes includes the subbands analog to 
monolayer or the bilayer graphene; 10 ' 13 and the Hamil- 
tonian is actually decomposed into independent subsys- 
tems effectively identical to monolayer or bilayer . 14 ' 16 
The ABC graphite has a quite different electronic struc- 
ture from ABA'si^^^i 9 .. In particular, the low- 
energy band of a finite ABC multilayer are given by the 
surface states localized at outer-most layers; 10 ' 15 and the 
interlayer potential asymmetry opens an energy gap in 
those bands i 26 ' 46 ' 49 This is in sharp contrast with ABA 
multilayers where potential asymmetry causes a band 
overlapping ! 18 ' 26 

In considering the interlayer potential asymmetry in- 
duced by an external electric field, it is essential to 
take into account screening effect, as done in bilayer 
graphenej 2 ^— and ABA multilayer s 18 ' 20 ' 50 . Experimen- 
tally the interlayer screening effect in the gate electric 
field was probed in thin graphite films ££ir£& Here we 
calculate the self-consistent band structure of ABA and 



ABC multilayers with various iV's in the presence of per- 
pendicular electric field. For ABA multilayers, we show 
that the electric field generally produces band overlap- 
ping, and the screening is shown to be linear to the 
field amplitude. In ABC multilayers, on the other hand, 
the low-energy surface band causes a strong non-linear 
screening effect through opening an energy gap. The pa- 
per is organized as follows: we present the effective mass 
models for ABA and ABC multilayers in Sec. UH and 
compute the band structure including the self-consistent 
screening effect in Sec. Mil The conclusion is given in 
Sec. EH 



II. EFFECTIVE HAMILTONIAN AND BAND 
STRUCTURE 



A. ABA multilayers 

We first consider a multilayer graphene with ABA 
stacking, composed of N layers of a graphene layers. We 
label A and B on i-th layer as Aj and Bi . In ABA stack- 
ing, the sites B\, A%, -B3, A4 ■ ■ ■ are arranged along ver- 
tical columns normal to the layer plane, while the rest 
sites Ai, 2?2, A3, B4 ■ ■ ■ are above or below the center of 
hexagons in the neighboring layers, as shown in Fig.[T](a). 
The system is described by a k-p Hamiltonian based on 
three-dimensional (3D) graphite modeli 3 ^— As the sim- 
plest approximation, we include parameter 70 describing 
the nearest neighbor coupling within each layer, and 71 
for the coupling of the interlayer vertical bonds. The 
band parameters were experimentally estimated in the 
bulk ABA graphite, for example 42 as 70 = 3.16 eV and 
71 = 0.39 eV, which we will use in the following calcu- 
lations. The lattice constant of honeycomb lattice (dis- 
tance between nearest A atoms) is given by a = 0.246 
nm, and the inter-layer spacing d = 0.334 nm. 

The low-energy spectrum is given by the states in the 
vicinity of K and K' points in the Brillouin zone. Let 



\Aj) and \Bj) be the Bloch functions at the K point, 
corresponding to the A and B sublattices, respectively, 
of layer j. If the basis is taken as \Ai), \Bi); \A 2 ), \B 2 ); 
\An) ,\Bn) , the Hamiltonian around the K point is 
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given by 
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B. ABC multilayers 
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where Uj is the electrostatic potential at jth layer, and 
we defined p± — p x ± ip y with p = — ifiV. v is the band 
velocity of monolayer graphene given by v — \/3ajo/2h. 
The effective Hamiltonian for another valley, K' , is ob- 
tained by interchanging p + and p- ^ 

When Uj =0, Hamiltonian (JTJ can be decomposed into 
subsystems identical to bilayer or monolayer graphenes 
with a basis appropriately chosen^ The subsystems are 
labeled by an index m which ranges as 



The ABC multilayer have a different arrangement 
shown in Fig. [T] (b), where vertical bonds couple the 
pairs (Bj,Aj + i) for j = 1, 2, • • • , N — 1. We use the 
same notation 70 and 71 as in ABA graphite, for the 
nearest intralayer and interlayer coupling, respectively. 
Although the band parameters are not identical between 
ABA and ABC graphites, we refer to the values of 
ABA in the following numerical calculations, assuming 
that the corresponding coupling parameters have similar 
values. 44 Hamiltonian around the K point can be written 

10.44.45.48 
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with the same matrices defined in Eq. ©. When Uj=0, 
the eigenenergies are given by 
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The eigenenergies at Uj = are given-i^ for m = as with s = ± &n d ip n (n = 1,2,- ■■, N) being solutions of 
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where p ■■ 



p^, pi = ±, s — ±, and 

7T m7r 
' m = 2 ~ 2(iV+ 1)' 



M 2 ,(4) 
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The corresponding wavefunction is \tp) = ip(Ai)\Ai) 
ipiBijlBt) + ■■■ with 



(l>(Aj)\ _ „ (e ie V-V sin(N + l-j)<p n 



(11) 



m = only exists in odd-layer graphene and gives an 
energy band identical to monolayer graphene. Other m's 
are bilayer-type band structures where fi = — gives a pair 
of electron (s — +) and hole bands (s = — ) touching at 
zero energy, and (J, = + another pair repelled away by 
±271 cos n m . The dispersion around k = is approxi- 
mately quadratic with the effective massi 



7i 

TO = — 5 COS K m , 
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giving the density of states at zero energy, p m = 
g v g s m* /(2irh), with = 2 and g s = 2 arc valley (K, K') 
and spin degeneracies, respectively. 

The quantity K m corresponds to the wave number k z 
in the layer stacking direction (z-direction) via K m = 
A: z cij 10 ' 14 The wave function of subband m is indeed a 
standing wave in z-direction with wave number n m . The 
total density of states per layer, p = (1/N) J2 m Pm: ap- 
proximates in large N limit, 



where 9 = a,rcta,np y /p x and C is a normalization factor. 
In the bulk limit, ip n corresponds to the wavenumbcr 
along the layer stacking (z) direction. Actually, Eq. QlOp 
is obtained by imposing a condition that a standing wave 
in z-direction, composed by bulk wavefunctions, becomes 
zero at fictitious sites B and A^+i out of the system. 

Equation l|10p has N solutions of ip giving independent 
eigenstates. All of ip n are real when vp > jiN/(N + 1), 
while only one becomes complex when vp < "f±N/(N+l), 
which corresponds to the evanescent mode in the bulk. 
In vp <C 71, the complex branch approximates e ltp w 
—vp/ji, giving the dispersion 
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with s = ±. These are the only bands which appear at 
e = and dominate the low-energy physics. The corre- 
sponding wavefunction is 



7i 



P ~ 9v9s 2^Vv^ 
where ^ m is replaced with integration in k. 
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The wave amplitude becomes largest on the top and bot- 
tom layers and decays exponentially inside, and thus is 



3 



(a) ABA 



(b) ABC 



Ai, 



A 2 



A30 



B 3 



TO 

B21 




To 

A1 c * B1 



A 2 



A3' 



B 2 



B 3 




At 



A2 f B2 



A3 <^—f B3 



n 



FIG. 1: Atomic structures of multilayer graphenes with (a) 
ABA (Bernal) stacking and (b) ABC (rhombohedral) stacking 



III. SCREENING EFFECTS 
A. Self-consistent treatment of screening effect 

We compute the band structure of ABA or ABC mul- 
tilayer graphenes in presence of gate electric field taking 
account of the screening effect. We consider undoped 
free-standing multilayer graphenes with an external elec- 
tric field Fq applied to the perpendicular direction. This 
situation can be realized in an experimental set up with 
an external top and bottom gates electrodes which are 
held at the opposite gate voltages with respect to the 
grapheneJ^ 

The potential at each layer, Uj(j = 1, 2, • • • , N) should 
be determined self-consistently. If a set of Uj is given, we 
can compute the band structure using the Hamiltonian 
Eq. (UJ for ABA or Eq. © for ABC multilayers. Then 
we determine the Fermi energy so that the total density 
is equal to n to t (= in the present case), and calculate 
the electron density at each layer, nj(j = 1,2, •• • ,N), 
from the occupied eigenstates. For screening effect, we 
consider the multilayer as parallel plates with zero thick- 
ness and respective electron densities rij. The electric 
field between jth and (j + l)th layers is then given by 



regarded as a surface stated. The wave function is ex- 
actly localized at the sites Ax and Bn at p = 0, and as p 
increases, the decay length increases as — 1/ log(wp/7i) in 
units of interlayer spacing d. In Fig. |4j we plot the band 
structures of ABC graphenes with N = 2,3, 5, 10 and 20, 
where the results of [/,- = are indicated as black dotted 



curves. The surface states of Eq. (jl2[) are shown as a pair 
of electron and hole bands touching at e = 0, which be- 
come flatter as N increases. The bilayer graphene (AB) 
can be regarded as N = 2 of ABA family and also that of 
ABC family, and indeed, equally described either of Eqs. 
© or ©. 

When we consider the low-energy physics around zero 
energy, it is convenient to use the effective Hamiltonian 
reduced to the basis | gjv)^ ' 15 ' 47 In presence of Uj, 
it reads 
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This approximation is valid when vp/ji <C 1, i.e., the 
actual wave function, Eq. (p~5|) . is well localized to A\ or 
Bn. When we set the origin of potential as U\ + Un — 0, 
the eigenenergy is given by 



s s , p = Sv / 7l 2 (W7i) 2Ar + (Af//2) 2 , (15) 

where AU = U\ — Un and s = ±. The potential differ- 
ence AU opens an energy gap between the valence and 
conduction bands. 
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Here e is the permittivity of the interlayer spaces with- 
out the screening effect of 7r-band electrons, and we set 
e = 2 in the following calculations. Eq. (| 16[) immedi- 
ately gives a new set of the electrostatic potential Uj, 
which should be identical to the initial Uj. To find the 
self-consistent solution, we employ an iterative numerical 
approach, where we start with Uj = eFo[j— (N + l)/2] as 
initial values and iterate the process until Uj's converge. 



B. ABA multilayers 

In Fig. [5] (a), solid red curves show the self-consistent 
band structures of ABA multilayers with several N's, in 
presence of the external field eF$d = 0.2 7 i. The original 
band structures at Fq = are also shown as dotted black 
curves. In N > 3, we see that the lowest electron band is 
pulled down and the highest hole band is lifted up, mak- 
ing a band overlap around zero energy, as was previously 
recognized in the case of TV = 3 and 4 . 18 ' 26 The energy 
width of overlap becomes almost constant in A^IO. 

Figures [2f b) and [2c) show the corresponding poten- 
tial distribution Uj and electron density rij , respectively, 
at the same external field eF$d = O.271. In N > 10, 
we observe that the electric field (i.e., gradient in Uj) 
is screened within a few layers from the surface, leav- 
ing a triangular potential pocket at each end. The po- 
tential decay near the edge is almost identical between 
N = 10 and 20. The overlapping bands observed in 
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(a) Band structure (ABA, eFod/y = 0.2) 
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(b) Potential (ABA, eF d/y, = 0.2) (c) Electron density (ABA, eFod/y, = 0.2) 
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FIG. 2: (a) Self-consistent band structures of ABA (Bernal) multilayer graphenes with several layer number N's, at external 
field eFod/^/i = 0.2 (red, solid) and (black, dotted), (b) Potential distribution and (c) electron density of ABA multilayers 
with several N's at eFod/ji = 0.2. 
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FIG. 3: Potential difference between outermost layers as a 
function of external field, in ABA multilayer graphenes with 
several N's. 

Fig. [2] (a) are actually the bound states trapped at ci- 
ther of pockets; the states of the lowest electron and the 
highest hole bands are indeed localized at the potential 
minimum (left end) and maximum (right), respectively. 



Since Ep is zero, those bands are populated by electrons 
or holes, contributing to the most part of the screening 
field. A smooth decay observed in the electron density 
appears different from Rcf. 50 , which finds a charge os- 
cillation with every second layer. We presume that this 
is due to the contribution from the intraband excitation, 
which was dropped in numerical calculations for neutral 
systems in Ref<££. 

The typical screening length A s (penetrating depth 
of electric field) can be roughly estimated by Thomas- 
Fermi approximation^ In this treatment, the potential 
decay on the surface is expressed as U{z) cx e~ z l Xs with 
X s = (e 2 p3D/e)~ 1 ^ 2 , where p 3 D is the three dimensional 
density of states at the Fermi energy. For graphene, if 
we substitute p^u — p/d with p of Eq. ([7]), we obtain 50 

Using the parameters above, we get A s ~ 1.3d rs 0.43nm. 
In Fig. [2] (b) , we plot an exponential curve with decay 
length X s in Eq. (IT71) as a dotted curve to fit with the 
right half of the curve of N — 20, which shows a fairly 
nice agreement. The depth of potential depth, or \U(z — 
0)|, is roughly estimated as eFo\ s , which determines the 
order of the energy width in band overlapping. 
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FIG. 4: (a) Self-consistent band structures of ABC (rhombohedral) multilayer graphenes with several layer number N's, at 
external field eFod/ji = 0.2 (red, solid) and (black, dotted), (b) Potential distribution and (c) electron density of ABC 
multilayers with several N's at eFod/yi — 0.2. 



Figure [3] displays the potential difference AU = U\ — 
Un as a function of the external field i*o- AU rises almost 
linearly in increasing Fq, except for a slight sub-linear 
components in large Fq . This is consistent with Thomas- 
Fermi approximation, since it gives linear screening in a 
weak external field. 



and bottom layers, 5n = n(A\) — h(-Bjv)> is calculated as 
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where (^_ iP (Ai),^_ jP (.Biv)) is the eigenvector of Eq. (JTH 
for s = — band, and 



C. ABC multilayers 



The screening property of ABC multilayers is quite 
different from that of ABA, as the density of states di- 
verges at e = due to the fiat band of the surface states. 
Before numerical calculations with full band model, we 
present an analytical approach using the effective 2x2 
Hamiltonian of Eq. (|14p valid in low energies. The po- 
tential difference AU between the top and bottom layers 
opens an energy gap between the valence and conduction 
bands, and thus only the lower band (s = — ) is occupied 
when ntot = 0. The density difference between the top 
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with r(a;) is the gamma function. 

The density imbalance Sn causes the screening field 
F inc j = —e5n/(2e) opposed to the external field Fq, re- 
sulting in the total potential difference AU = e(Fo + 
Fi a d)(N - l)d. Together with Eq. ([15]). we obtain the 
self-consistent equation for AU, 



AU=e(N-l)d 
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In N > 3, AU is negligible compared to AU 2 ^ N when 
AU is small enough. Then the equation is solved ap- 
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FIG. 5: Potential difference between outermost layers of ABC 
multilayers as a function of external field. Vertical broken line 
indicates the critical field F c . Lower four panels compare the 
same results (solid curves) to the approximate expressions of 
Eq. (PU (dashed). 



proximately as 



AU w 2 7lJ Fo 
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which is essentially non- linear in Fq. In large- N limit, 
we have f N &l/2 and thus AU w 27i(F /.F c ) JV / 2 , where 
F c = en c / (2e) is a characteristic field with an associated 
electron density 
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In increasing Fq, AU rapidly increases when the external 
field exceeds F c . 

n c is the electron density accommodated in the flat- 
band region in large N limit (up/71 < 1), i.e., the number 
of surface states. The field is completely screened in F n < 
F c , because the surface states are able to supply positive 
and negative charge to opposite surfaces to cancel the 
external field. The screening collapses at F c , when the 



density required for canceling exceeds the surface states 
population n c . 

N = 2 (AB) is an exceptional in that the integration 
in Eq. (|19p diverges logarithmically, giving infinite 5n. 
Actually this is an artifact of the reduced 2x2 model, 
due to the incorrect contributions from large p where the 
reduced Hamiltonian is not accurate. We can remove this 
by introducing a momentum cut-off p c ~ Ji/v, and get 
/n—2 ~ (1/2) log(7i/A[/). When we neglect the loga- 
rithmic dependence of f^, AU becomes linear in Fq in 
accordance with Eq. (f2"U| . The logarithmic factor gives a 
weak singularity at AU = 0. 

Now we numerically calculate the self-consistent band 
structure of ABC multilayers using the full Hamiltonian 
Eq. (J5]). Figure S] (a) shows the results at eF d — O.271 
(red, solid) and (black, dotted). In presence of the 
external field, an energy gap opens at low-energy as ex- 
pected. The gap width becomes smaller in N > 5 in in- 
creasing N, suggesting a strong screening effect in large 
stacks. Figures IUb) and [He) show the corresponding 
potential distribution Uj and the electron density rij , re- 
spectively, at the same field eF d = 0.2 7 i. At N = 20, 
the potential is almost flat inside, as the external field is 
mostly screened by the electric charge on surface states 
localized at the outermost layers. This is in contrast with 
ABA multilayers, where an external field always pene- 
trates inside with a few-layer thickness. 

Figure [5] shows the plots of the potential difference 
AU = U\ — Un as a function of the external field Fo. We 
actually observe non-linear behavior expected in the ana- 
lytical argument, where AU rapidly increases at F ~ F c 
(shown as a dashed vertical line). Lower panels in Fig. 
[5] compare the numerical results (solid) to the analytical 
expression Eq. (|2Tj) (dashed). We have nice agreements 
for N < 5 in small Fq, while the approximation becomes 
worse in large stack of N > 10. In large N's, the low- 
energy band almost reaches vp/ji ~ 1, where the wave 
function deeply penetrates into the bulk in accordance 
with Eq. (jT5|) . The finite penetration length makes the 
screening less effective, compared to the previous 2x2 
model assuming the wave functions perfectly localized 
on the surface layers. As a result, numerical curves at 
large iV's rise less sharply than the analytical ones as ob- 
served in Fig. O The wave penetration to the bulk is 
also responsible for Mexican hat structure^ or narrow- 
ing of the gap around vp/ji ~ 1 observed in Fig. [4] (a). 
There the actual energy splitting becomes smaller than 
AU because the wave function is not perfectly localized 
at surface layers. 

The width of the energy gap is an important quantity 
which can be detected experimentally. Figure[H]shows the 
gap width against the external field in the self-consistent 
band structures of ABC multilayers. In N < 5, the band 
bottom is approximately flat and the gap width therefore 
approximates AU (the splitting at p — 0), and actually 

rises in proportional to F^ 2 . In large stacks of N > 10, 
the energy gap becomes maximum around Fq ~ F Cl and 
is suppressed in greater -Fb's, due to the gap narrowing 
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FIG. 6: Energy gap width in the self-consistent band struc- 
tures as function of external field, for ABC multilayers with 
several iV's. Vertical broken line indicates the critical field 
F c . 



ABA multilayers, the electric field produces band over- 
lapping accompanying a linear screening well described 
by the Thomas-Fermi approximation. In ABC multi- 
layers, in contrast, the surface state bands dominating 
low energies cause a strong non-linear screening effect 
through opening an energy gap. 

While in the present model we only include the primary 
parameters 70 and 71 in our model, the extra band pa- 
rameters corresponding to the further hopping generally 
affect the band structure of multilayer graphenesi 36 ' 38 " — 
In ABC graphenes, it was shown that the extra param- 
eters gives a fine structure to the surface band, of which 
energy scale is expected to be of the order of 10 meV.— 
We expect that the screening property would be influ- 
enced by those effects when the external potential is as 
small as those energy scales. As another remark, the 
electron-electron interaction other than the screening ef- 
fect may create non-trivial ground states in a flat band 
such as in ABC multilayers, while we leave those prob- 
lems for future works. 



at vp/71 ~ 1. 

IV. CONCLUSION 

We studied electronic band structures of ABA and 
ABC graphcnc multilayers in the presence of an perpen- 
dicular electric field, including the screening effect. In 
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